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Two sorts of bosons in an optical lattice at commensurate filling factors can form five stable 
superfluid and insulating groundstates with rich and non-trivial phase diagram. The structure of 
the groundstate diagram is established by mapping d-dimensional quantum system onto a (d + 1)- 
dimensional classical loop-current model and Monte Carlo (MC) simulations of the latter. Surpris- 
ingly, the quantum phase diagram features, besides second-order lines, first-order transitions and 
two multi-critical points. We explain why first-order transitions are generic for models with pair- 
ing interactions using microscopic and mean-field (MF) arguments. In some cases, the MC results 
strongly deviate from the MF predictions. 
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Ultracold atoms trapped in an optical lattice (OL) 
form an intriguing strongly correlated quantum 
system. The unprecedented control over parameters of 
the effective Hubbard-type Hamiltonian renders this sys- 
tem an important object for the study of quantum phase 
transitions 0. Single-component bosons without inter- 
nal degrees of freedom have only two phases in a regu- 
lar lattice: superfluid (SF) and Mott-insulator (MI) (at 
commensurate filling factor |^). When several bosonic 
species are combined in the OL, the naive expectation 
that their groundstates are straightforward mixtures of 
MI and SF with respect to participating components is 
wrong — the phases of spinor and multi-component sys- 
tems are far more subtle 000001111. 

In this Letter, we study a commensurate two- 
component bosonic system described by the on-site Hub- 
bard Hamiltonian: 

H = - ^ {t a a\ a a, ja + H.c.) + i U aa - n ia n ial . (1) 



<ij> cj 



Here at creates a boson of the sort a = A,B on site 

i, riic — a\ a a ia , and < ij > denotes pairs of nearest- 
neighbor sites. In what follows, we consider only equal 
filling factors of the components, tia = ub = n, with n 
integer, and, for brevity, denote U ab = — V, U aa = U a . 
Similar (incommensurate) two-species bosonic model has 
been studied recently to look at the differences with the 
fermionic Hubbard model 

At double commensurate filling, recent mean-field 
analysis of the model in ^2] failed to reveal phases 
and, correspondingly, phase transitions which can not be 
reduced to simple mixtures of single-component states. 
This conclusion is disappointing considering predictions 
of other strongly correlated superfluid groundstates for 
two-component inconvertible bosons: a paired super- 
fluid vacuum (PSF), which is equivalent to the superfluid 
state of diatomic molecules and to BCS superconductor 
000.01; and a super-counter-fluid (SCF), in which 



the net atomic superfluid current is zero, and yet the 
equal currents of the components in opposite directions 
arc superfluid 013. 

In this Letter, we perform MC simulations of the 
(d + I)-dimensional classical analog of the on-site Hub- 
bard model Q and find five stable superfluid and insu- 
lating phases: (i) MI, (ii) MI of sort A and SF of sort B 
(MIb + SF A ) and its A <-> B analog, (iii) SF of sort A and 
SF of sort B (2SF), (iv) PSF, and (v) SCF. An interact- 
ing mixture of two mutually penetrable superfluids (2SF) 
exists even without the optical lattice and corresponds to 
the t ^> U a<(T i limit. This state has two non-zero complex 
order-parameters (tPa) and {ips)- By increasing cither 
Ub or V a one drives the corresponding component from 
the superfluid to the Mott-insulating state; accordingly, 
in MIb + SFa we have finite (iPa), and zero (tps)- When 
both U a and Ub are strong, the groundstate is MI with 
all order parameters being zero. Our proof then concerns 
the existence of PSF and SCF phases; in both phases 
(i/>a) = {^b) = 0, while <t>PSF = (Mb) ± in PSF, 
and &scf = (ipAips) ^ in SCF. It is worth noting 
that, while the PSF, representing atomic A+B pairing, 
requires V > 0, the SCF describes pairing of particles A 
and holes B and occurs when V < 0. 

The most surprising MC result is that 2SF-MI tran- 
sition may be I-order. This result also follows from the 
mean field (MF) analysis of the problem along the lines 
suggested in for the single-component case. Finally, 
we develop microscopic arguments explaining why the I- 
order transition is generic for models with pairing inter- 
actions, and show that MI groundstates may be further 
classified in terms of their excitation spectrum. 

To prove that possible groundstates of Eq. include 
PSF, we assume the inter-exchange symmetry A < — > B, 
implying t& = tg = t and Ua — Ub = U , and consider 
the limit described by two strong inequalities: t/U <C I 
and 7/Z7 <C I, with 7 = U - V [here V > and 7 > 0; 
at 7 < the system collapses.] Then, the effective low- 
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energy Hilbert space is determined by states were on each 
site HiA = niB] these are separated from other states by 
a large pair-breaking gap « U. We thus naturally arrive 
at the description in terms of pairs. In the second-order 
perturbation theory in t/U (cf., e.g., @) the dynamics 
of pairs is given by the effective Hamiltonian (we omit 
terms proportional to the total number of particles): 

H v = -t \{°i°t + H.c.) + 2m i m j ] + 7 ^ m\ . (2) 

<ij> i 

Here rrii are pair occupation numbers, the raising oper- 
ator Of is defined by (m'^Oflmi) = (rrii + l)<5 m ' imi +i, 
0~[ = (Of y, and t — 2t 2 /U. In contrast to the standard 
single-boson hopping that scales linearly with the typi- 
cal occupation number the hopping amplitude for pairs is 
quadratic in rrii ■ If potential energy terms in Eq. J2J) were 
omitted, the groundstate would collapse to a droplet with 
the diameter comparable to the lattice constant. The 
second term in the brackets, Eq. (J2J), describes nearest- 
neighbor attraction, and further enhances collapse insta- 
bility. A stable groundstate arises only when the on-site 
repulsion 7 is strong enough. On the other hand, at very 
large 7 the commensurate groundstate is MI. The ques- 
tion is then whether for some t/7, the groundstate is PSF 
rather then MI or collapsed. The positive answer is read- 
ily seen in the limit of very large molecular filling factor 
m = n/2 ^> 1. In the region i«7< m 2 t the system is 
stable against collapse; the nearest-neighbor attraction is 
negligible. Since the maximum insulating gap ~ 7 does 
not depend on m 0, we conclude that for 7 < m 2 t the 
groundstate must be superfluid. [Using mapping to the 
quantum rotor model [3] we know, in fact, that MI state 
requires 7 > m 2 t\. Finally, for m — 1 we have performed 
quantum MC simulations of Eq. in d = 2 and found 
that at 7 = lOt the groundstate is superfluid. 

Similarly, the existence of the SCF phase can be shown 
in the limit t/U < 1 and \j\/U < 1 [where 7 = U + V 
and V is negative] studied in Ref. Q. The effective 
Hamiltonian can be now written in terms of the spin- 
s' = (n A + n B )/2 operators, H s = -t J2<ij> S » S j + 
1^2i{Siz) 2 , which for small positives has the easy-plane 
ferromagnetic groundstate, or SCF p|. Furthermore, one 
can show that PSF and SCF are qualitatively similar and 
SCF may be viewed as the result of pairing between parti- 
cles of one component and holes of another |l3j . Because 
of this equivalence, we discuss below only the V > case, 
that is, the case of the PSF. 

To reveal the global structure of the phase diagram 
we performed Monte Carlo simulations for the (d + 1)- 
dimensional classical analog of the bosonic Hubbard 
model. The so-called J-current model [l4| is built on 
particle worldlines (space-time currents) in discrete imag- 
inary time, and we straightforwardly generalize it to the 
two-component case: 

a, a' i 



Here are integer-value currents [(d+ l)-dimensional 
vectors] subject to the local zero-divergence constraint, 

V- ■// — 0, and U aa > ~ U aa > jt relate the effective action 
parameters to the original Hubbard Hamiltonian. This 
model has the same superfluid and insulating phases as 
Eq. QJ, and we use it to understand the topology of 
phase boundaries, the existence of multicritical points, 
and first-order lines. We find it convenient to fix U a and 
Ub and to plot results in the (V, f)-plane, where v ~ 1/t 
is the scaling factor for all three dimensionless param- 
eters. The superfluid phases are identified by looking 
at various superfluid stiffnesses, pf^ = ([W ( ^] 2 ) /dL d ~ 2 , 
p (PSF/scF) = ± wW] 2 )/DL D - 2 , expressed in 

terms of the winding number fluctuations [15|, where 
tfr(v) = L -i ^ j} a) (superfluid stiffness and compress- 
ibility are equal in the space-time symmetric model). 

In Fig. ^ we present the phase diagram of the two- 
component J-current model in (d + 1 = 3)-dimensions. 
Corresponding superfluid stiffness goes to zero continu- 
ously when approaching the lines of the critical points 
labeled as U(l). The correlation radius exponent (ob- 
tained from finite-size corrections) is consistent with the 
known value for the U(l) universality class in 3D. The 
first-order transition was identified by (i) double-peak 
structure of the energy distribution function (in small- 
size systems), and (ii) hysteresis loops in all quantities 
(in large-size systems). Though we have not performed 
exhaustive MC study of the phase diagram in other di- 
mensions we found (i) the I-order 2SF-MI transition in 
d = 3, and (ii) no evidence for the first-order 2SF-MI 
transition in d = 1. 

The I-order 2SF-MI line in the symmetric case ( Ua — 
Ub), becomes strongly suppressed by the anisotropy, 
V a~Ub 7^ 0, between the components. For Ub/Ua = 2, 
the point where all four phases meet is already a simple 
cross of two U(l) lines — decoupled U(l) x U(l) tetra- 
critical point (see, e.g., 0])- Points where the I-order 
line starts and ends represent multicritical points. 

Normally, first-order transitions can be qualitatively 
accounted for in simple mean-field models, and we pro- 
pose such a model for our case. Away from the multi- 
critical region, all transitions are of the [/(l)-universality 
class and, thus, described by the corresponding l^l 4 ac- 
tions 3 for atomic, tpA> V'fl) an d molecular, fields. We 
arrive at the simplest effective free energy by combining 
three l^l 4 actions and writing the interaction term in the 
form of the molecule "creation/annihilation" process out 
of A- and B-particles. Omitting gradient terms: 

^ = \[rA\4>A\ 2 +r B \4> B \ 2 + r M \<5>\ 2 ] +1[ 5A |^| 4 

+ 9b\^b\ 4 +. 9M |$| 4 ] -g(**MB+cc) ■ (4) 

The mean-field phase diagram follows from minimiza- 
tion of T . In Fig|3 taking advantage of the scaling free- 
dom for all the fields and T , wc set g a = §b = 9m — 
g = 1. It reproduces correctly the topology of bound- 
aries between the phases and some of their properties. If 
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FIG. 1: Phase diagrams of the d = 2 two-component J-current model in the (V, i^)-plain for the symmetric, Ua = i/s = 2 (left), 
slightly asymmetric, 17a = 1-6, Ub = 2 (center), and strongly asymmetric, Ua = 1, t/s = 2 (right) models. The l st -order phase 
transition line is dramatically reduced in the presence of weak anisotropy, and completely disappears for strong anisotropy 
between the components. All horizontal errorbars are smaller than point sizes (typically of order 10 -3 ), and lines are used to 
guide the eye and to distinguish between different phase boundaries. The insert shows more clearly the region where 2SF, PSF 
and MI phases meet. (Commensurability and intrinsic symmetry of the J-current model result in a straightforward mapping 
of the SCF regime onto PSF one: V -» - V.) 
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FIG. 2: The mean-field diagram for va = r, rs = r + 0.5 The 
bold solid line corresponds to the first-order transition, and 
dashed lines describe continuous transitions. In the limit of 
tb — ► va, the SFa + MIb domain vanishes. 



the A-B asymmetry is not large, it features a I-order line, 
see Fig. |21 In the strongly anisotropic case, |rg — ta\ > 1, 
the MF theory also captures the disappearance of the I- 
order transition. On another hand, the prediction of the 
pronounced I-order 2SF-PSF line does not agree with 
the numerical data in 3D. With our system sizes up to 
128 3 sites we did not find any evidence for the I-order 
2SF-PSF or 2SF-(MI B + SF A ) transitions, see insert in 
Fig. ^ It is probably too early to draw the final conclu- 
sion on the structure of the multicritical point, because 
a similar study in 4D (for the system size 32 4 ) revealed 
a tiny (but finite) I-order 2SF-PSF line. In any case, the 
suppression of the I-order 2SF-PSF transition constitutes 
a strong deviation from the MF prediction. 

First-order SF-MI transition in the single- component 
system. It is generally accepted that in the single- 
component, commensurate Bose system the SF-MI tran- 
sition is continuous |J] . Numerous simulations of the on- 
site Hubbard and J-current models perfectly agree with 
this picture (for the latest simulation see [13]). 



Excitations in MI are gapped and described as quasi- 
particles and quasiholes with the relativistic disper- 
sion law at small momenta (for small gaps): e(k) = 
\/A 2 + c 2 fc 2 , where c is the velocity of sound in the 
SF phase. The dilute gas of quasiparticles is charac- 
terized by the effective mass m* = A/c and some s- 
wave scattering amplitude, a* (to be specific, we as- 
sume that d = 3). If the scattering length is positive, 
the state of the dilute excitation gas with density n qp 
is stable, and the energy density cost of creating it is 
E qp = (A — fi)riqp + (27ra*/m, t )n 2 lp , where fi is the chem- 
ical potential. Since the effective longwave action for the 
U(l)-transition has positive coefficient in front of the |^| 4 
term, in the vicinity of the critical point the MI state is 
always described by positive a* . As the chemical poten- 
tial is increased above the threshold value A, the system 
state becomes superfiuid (this continuous MI-SF transi- 
tion, induced by adding extra particles/holes, is mean- 
field like 0). 

Imagine now a MI state with gaped quasiparticle 
excitations, but now with negative effective scattering 
length. Although the MI vacuum itself may remain sta- 
ble, the state of the quasiparticle gas at any small den- 
sity n qp |a*| 3 <C 1 is unstable against collapse to a dense 
droplet. We thus conclude that this MI will undergo a 
first-order phase transition to the superfiuid state at some 
/i. = |U C finite distance below A to gain negative poten- 
tial energy. Furthermore, if for some system parameters 
A ^ 0, but fi c = 0, the MI-SF transition in the commen- 
surate system will happen by I-order scenario too. 

We now argue that MI states with negative a* nat- 
urally arise in models with strong pairing interactions 
when potential energy favors two bosons on the same 
site, but increases fast for occupation numbers > 2 to 
prevent collapse. For sufficiently strong pairing, one may 
have a superfiuid state of tight molecules, or PSF (cf. 
When repulsion between molecules is increased, 
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FIG. 3: Sketch of possible MI phases in the commensurate, 
n =even, single-component Hubbard model with pairing in- 
teractions and phase transition lines from MI to the SF and 
PSF phases. The region between crosses is characterized by 
first-order MI-SF transition as a function of /i. 

PSF undergoes a standard continuous PSF-MI transi- 
tion. In the vicinity of the transition point the lowest 
excitations above MI are bound bosonic pairs. We ob- 
serve then, that depending on the value of the pairing 
interaction there should exist at least three MI ground- 
states distinguished by the value of the effective scatter- 
ing length of its quasiparticles: as we go along the line C 
in Fig. a* starts from positive value (the quasiparticlc 
gas is stable), then changes sign and becomes negative 
(the quasiparticle gas is unstable against collapse), and 
finally goes through the pole and changes sign again (the 
quasiparticle gas of molecules is stable). It seems unlikely 
that SF, PSF, and three different MI phases meet at the 
same point. In our view, the intersection of the SF-MI 
and | a* | = lines marks the beginning of the first-order 
SF-MI transition, while the intersection of the PSF-MI 
and |a*| = oo lines marks its end. In this picture, the 
critical point where the SF-MI line changes from con- 
tinuous to first-order is characterized by the continuous 
Lorentz-invariant action with zero term. 



In d = 2, the weak logarithmic dependence of a* on 
quasiparticle density does not change the qualitative pic- 
ture, because the first-order transition involves finite par- 
ticle density jumps. In d = 1, the notion of the scattering 
length is ill-defined, and two quasiparticles in the long- 
wave limit either form a bound state or repel each other 
like hard-core spheres. We conjecture then, that in d = 1 
(i) MI with first-order transition in fi does not exist, (ii) 
the SF-MI transition is always continuous. 

The above considerations readily generalize to the A- 
B symmetric two-component case. Now, the criterion 
for the MI groundstate, which is unstable against light 
doping by A- and B-particles (for the symmetric case 
A a = Ab = A), follows from the scattering matrix 
(a*)<T<T' becoming non-positive definite due to increasing 
attraction between the components. We note, that this 
criterion must be satisfied in a finite region in param- 
eter space since existence of AB-molecules implies that 
(o-*)ab can be arbitrarily large and negative before going 
through the pole corresponding to the formation of the 
bound state. This consideration, in complete analogy 
with the single-component case, suggests that PSF-MI 
and 2SF-MI lines are "bridged" by the first-order line in 
agreement with MC simulations and MF analysis. 

To explain the suppression and disappearance of the 
first-order region when the symmetry between the A and 
B components is broken (see Fig. we suggest that, 
for strong anisotropy, the lowest excitations above the 
MI groundstate in the whole parameter range are either 
A-particles or AB-molecules, and there is no reason for 
collective instability in the quasiparticle gas. Formally, 
this corresponds to pushing the |(a*)As| = oo line into 
the SFa + MIb phase — this possibility does not exist in 
the symmetric case. 
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